Studies 1-4: Figures

Luhrmann, Weisman, et al.

In this notebook, we create figures for Studies 1-4.

source("./scripts_general/dependencies.R")
source("./scripts_general/custom_funs.R")
source("./scripts_general/var_recode_contrast.R")
source("./study1/scripts_s1/s1_var_groups.R")
source("./study2/scripts_s2/s2_var_groups.R")
source("./study3/scripts_s3/s3_var_groups.R")
source("./study4/scripts_s4/s4_var_groups.R")

setwd("./study1/analysis/")
source("../../scripts_general/data_load.R")

d_all <- d1 %>%
  select(study, country, site, religion, subject_id, 
         pv_score, abs_score, spev_score) %>%
  mutate(religion = recode_factor(religion,
                                  "indigenous" = "Indigenous Religion",
                                  "charismatic" = "Charismatic Christianity"),
         # rescale to 0-1
         pv_score = pv_score/3) %>%
  full_join(d2 %>% 
              select(study, country, subj, 
                     abs_score, spev_score, dse_score) %>% 
              rename(subject_id = subj) %>%
              # rescale to 0-1
              mutate(spev_score = spev_score/4,
                     dse_score = dse_score/5,
                     religion = "General Population")) %>%
  full_join(d3 %>% 
              select(study, epi_ctry, epi_sample, epi_subj, 
                     por_score, spirit_score) %>%
              rename(country = epi_ctry, 
                     religion = epi_sample,
                     subject_id = epi_subj,
                     spev_score = spirit_score) %>%
              mutate(religion = recode_factor(religion,
                                              "general population" = "General Population",
                                              "charismatic" = "Charismatic Christianity"))) %>%
  full_join(d4 %>%
              select(study, p7_ctry, p7_subj, 
                     por_score, pv_score, abs_score, spev_score, dse_score) %>%
              rename(country = p7_ctry,
                     subject_id = p7_subj) %>%
              # rescale to 0-1
              mutate(por_score = por_score/2,
                     pv_score = pv_score/3,
                     spev_score = spev_score/4,
                     dse_score = dse_score/5,
                     religion = "General Population")) %>%
  mutate(religion = factor(religion, 
                           levels = c("General Population", 
                                      "Indigenous Religion",
                                      "Charismatic Christianity")),
         study = gsub("study", "Study", study)) 

Summary statistics

By study

Study 4

rsq_fun <- function(reg){
  reg_class <- class(reg)
  res <- rsquared(reg)
  
  if (grepl("lmerMod", reg_class)) {
    rsq <- res$Conditional
  } else {
    rsq <- res$R.squared
  }
  return(rsq)
}

d_sum_s4 <- d_all %>%
  filter(study == "Study 4") %>%
  group_by(country) %>%
  summarise_at(vars(spev_score, dse_score, pv_score, por_score, abs_score),
               funs(mean = mean(., na.rm = T), sd = sd(., na.rm = T))) %>%
  ungroup()

d_rsq_s4 <- data.frame(
  var = c("spev_score", "dse_score", "pv_score", "por_score", "abs_score"),
  rsq_fix = c(rsq_fun(lm(spev_score ~ country, d_all %>% filter(study == "Study 4"))),
              rsq_fun(lm(dse_score ~ country, d_all %>% filter(study == "Study 4"))),
              rsq_fun(lm(pv_score ~ country, d_all %>% filter(study == "Study 4"))),
              rsq_fun(lm(por_score ~ country, d_all %>% filter(study == "Study 4"))),
              rsq_fun(lm(abs_score ~ country, d_all %>% filter(study == "Study 4")))),
  rsq_ran = c(rsq_fun(lmer(spev_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 4"))),
              rsq_fun(lmer(dse_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 4"))),
              rsq_fun(lmer(pv_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 4"))),
              rsq_fun(lmer(por_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 4"))),
              rsq_fun(lmer(abs_score ~ 1 + (1 | country), 
                         d_all %>% filter(study == "Study 4")))))
d_rsq_s4

plot_s4_spev <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  geom_text(data = d_sum_s4 %>%
              mutate_at(vars(-country), funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(spev_score_mean, "\n(", spev_score_sd, ")")),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Spiritual Events")

plot_s4_dse <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = dse_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(dse_score, na.rm = T),
                              sd = sd(dse_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  geom_text(data = d_sum_s4 %>%
              mutate_at(vars(-country), funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(dse_score_mean, "\n(", dse_score_sd, ")")),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Daily Spiritual Experiences")

plot_s4_pv <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = pv_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(pv_score, na.rm = T),
                              sd = sd(pv_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  geom_text(data = d_sum_s4 %>%
              mutate_at(vars(-country), funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(pv_score_mean, "\n(", pv_score_sd, ")")),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Porosity Vignettes")

plot_s4_por <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = por_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(por_score, na.rm = T),
                              sd = sd(por_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  geom_text(data = d_sum_s4 %>%
              mutate_at(vars(-country), funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(por_score_mean, "\n(", por_score_sd, ")")),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Porosity Scale")

plot_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = country, y = abs_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(abs_score, na.rm = T),
                              sd = sd(abs_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  geom_text(data = d_sum_s4 %>%
              mutate_at(vars(-country), funs(format(round(., 2), nsmall = 2))),
            aes(y = 1, label = paste0(abs_score_mean, "\n(", abs_score_sd, ")")),
            color = "black", size = 2.5, vjust = 1) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Absorption")

fig_2_title <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 35))

fig_2 <- plot_grid(
  fig_2_title,
  plot_grid(plot_s4_spev, plot_s4_dse, 
            plot_s4_pv, plot_s4_por, 
            plot_s4_abs, NULL, 
            ncol = 2, labels = c("A", "B", "C", "D", "E")),
  ncol = 1, rel_heights = c(1, 20))

Removed 2 rows containing missing values (geom_point).

# fig_2

fig_2

Spiritual Events scores

Study 1

d1 %>%
  ggplot(aes(x = country, y = spirit_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spirit_score, na.rm = T),
                              sd = sd(spirit_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
  labs(x = "Country", y = "Study 1: Spiritual Events score (range: 0-1)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 2

d2 %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 2: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 3

d3 %>%
  # correct for scaling in original dataset
  mutate(spirit_score = spirit_score * 4) %>%
  ggplot(aes(x = epi_ctry, y = spirit_score, color = epi_ctry, fill = epi_ctry)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(epi_ctry) %>%
                    summarise(mean = mean(spirit_score, na.rm = T),
                              sd = sd(spirit_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 3: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

Study 4

d4 %>%
  ggplot(aes(x = p7_ctry, y = spev_score, color = p7_ctry, fill = p7_ctry)) +
  geom_jitter(height = 0, width = 0.25, alpha = 0.2, show.legend = F ) +
  geom_pointrange(data = . %>%
                    group_by(p7_ctry) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd),
                  shape = 23, color = "black",
                  show.legend = F) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  scale_y_continuous(limits = c(NA, 4), breaks = seq(0, 4, 1)) +
  labs(x = "Country", y = "Study 4: Spiritual Events score (range: 0-4)",
       caption = "Error bars are ±1 standard deviation from the mean")

All studies

d_all %>%
  ggplot(aes(x = country, y = spev_score, color = country, fill = country,
             group = study)) +
  geom_point(position = position_jitterdodge(jitter.width = 0.25,
                                             jitter.height = 0,
                                             dodge.width = 0.75), 
             alpha = 0.1) +
  geom_pointrange(data = d_all %>%
                    group_by(study, country) %>%
                    summarise(mean = mean(spev_score, na.rm = T),
                              sd = sd(spev_score, na.rm = T)) %>%
                    ungroup(),
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, shape = study), 
                  position = position_dodge(width = 0.75),
                  color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom") +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       shape = "Study", 
       caption = "Error bars are ±1 standard deviation from the mean")

d_spev_sum <- d_all %>%
  filter(!is.na(spev_score)) %>%
  group_by(study, country, religion) %>%
  summarise(mean = mean(spev_score, na.rm = T),
            sd = sd(spev_score, na.rm = T),
            n = n()) %>%
  ungroup()

Factor `religion` contains implicit NA, consider using `forcats::fct_explicit_na`

d_all %>%
  ggplot(aes(x = country, y = spev_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(~ study, scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_spev_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_spev_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Spiritual Events score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

Porosity scores

d_por_sum <- d_all %>%
  gather(por_scale, por_score, c(pv_score, por_score)) %>%
  mutate(por_scale = recode(por_scale,
                            "pv_score" = "Porosity Vignettes",
                            "por_score" = "Porosity Scale")) %>%
  filter(!is.na(por_score)) %>%
  group_by(study, country, religion, por_scale) %>%
  summarise(mean = mean(por_score, na.rm = T),
            sd = sd(por_score, na.rm = T),
            n = n()) %>%
  ungroup()

Factor `religion` contains implicit NA, consider using `forcats::fct_explicit_na`

d_all %>% 
  gather(por_scale, por_score, c(pv_score, por_score)) %>%
  mutate(por_scale = recode(por_scale,
                            "pv_score" = "Porosity Vignettes",
                            "por_score" = "Porosity Scale")) %>%
  filter(!is.na(por_score)) %>%
  ggplot(aes(x = country, y = por_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(cols = vars(study, por_scale), scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_por_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_por_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Porosity score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

Absorption scores

d_abs_sum <- d_all %>%
  filter(!is.na(abs_score)) %>%
  group_by(study, country, religion) %>%
  summarise(mean = mean(abs_score, na.rm = T),
            sd = sd(abs_score, na.rm = T),
            n = n()) %>%
  ungroup()

d_all %>% 
  filter(!is.na(abs_score)) %>%
  ggplot(aes(x = country, y = abs_score, 
             color = country, fill = country,
             group = religion)) +
  facet_grid(. ~ study, scales = "free", space = "free") +
  geom_point(position = position_jitterdodge(jitter.width = 0.8,
                                             jitter.height = 0.02,
                                             dodge.width = 0.75), 
             alpha = 0.15) +
  geom_pointrange(data = d_abs_sum,
                  aes(y = mean, ymin = mean - sd, ymax = mean + sd, 
                      shape = religion), 
                  position = position_dodge(width = 0.75),
                  fill = "black",
                  color = "black") +
  geom_text(data = d_abs_sum %>%
              mutate(ypos = case_when(
                grepl("charismatic", tolower(religion)) ~ mean + sd + 0.05,
                TRUE ~ mean - sd - 0.05)),
            aes(y = ypos, label = paste0("n=", n)), 
            position = position_dodge(width = 0.75),
            size = 3, color = "black") +
  scale_shape_manual(values = 21:24) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  guides(color = F, fill = F, 
         shape = guide_legend(override.aes = list(fill = "black"))) +
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "Country", y = "Absorption score (rescaled to 0-1)",
       # caption = "Error bars are ±1 standard deviation from the mean",
       shape = "Religion")

r1_spev <- lm(spirit_score_std ~ country, d1)
r1_pv <- lm(por_score_std ~ country, d1)
r1_abs <- lm(abs_score_std ~ country, d1)

r2_spev <- lm(spev_score_std ~ country, d2)
r2_dse <- lm(dse_score_std ~ country, d2)
r2_abs <- lm(abs_score_std ~ country, d2)

r3_spev <- lm(spirit_score_std ~ epi_ctry, d3)
r3_por <- lm(por_score_std ~ epi_ctry, d3)

r4_spev <- lm(spev_score_std ~ p7_ctry, d4)
r4_dse <- lm(dse_score_std ~ p7_ctry, d4)
r4_por <- lm(por_score_std ~ p7_ctry, d4)
r4_pv <- lm(pv_score_std ~ p7_ctry, d4)
r4_abs <- lm(abs_score_std ~ p7_ctry, d4)

df1 <- data.frame(study = "study 1",
                  var = c("spiritual experience", "porosity", "absorption"),
                  scale = c("spiritual events", "porosity vignettes", "absorption"),
                  rsq = c(rsquared(r1_spev)$R.squared, 
                          rsquared(r1_pv)$R.squared, 
                          rsquared(r1_abs)$R.squared))

df2 <- data.frame(study = "study 3",
                  var = c("spiritual experience", "spiritual experience", "absorption"),
                  scale = c("spiritual events", "DSE", "absorption"),
                  rsq = c(rsquared(r2_spev)$R.squared, 
                          rsquared(r2_dse)$R.squared, 
                          rsquared(r2_abs)$R.squared))

df3 <- data.frame(study = "study 2",
                  var = c("spiritual experience", "porosity"),
                  scale = c("spiritual events", "porosity scale"),
                  rsq = c(rsquared(r3_spev)$R.squared, 
                          rsquared(r3_por)$R.squared))

df4 <- data.frame(study = "study 4",
                  var = c("spiritual experience", "spiritual experience", 
                          "porosity", "porosity", "absorption"),
                  scale = c("spiritual events", "DSE", 
                            "porosity scale", "porosity vignettes", "absorption"),
                  rsq = c(rsquared(r4_spev)$R.squared, 
                          rsquared(r4_dse)$R.squared,
                          rsquared(r4_por)$R.squared, 
                          rsquared(r4_pv)$R.squared, 
                          rsquared(r4_abs)$R.squared))

df_all <- full_join(df1, df2) %>% full_join(df3) %>% full_join(df4) %>%
  mutate(var = factor(var, 
                      levels = c("spiritual experience", 
                                 "porosity", "absorption"))) %>%
  select(var, scale, study, rsq) %>%
  arrange(var, scale, study)

df_all %>% 
  mutate(percent_exp = paste0(round(rsq * 100), "%")) %>%
  select(-rsq) %>%
  spread(study, percent_exp) %>%
  mutate_at(vars(starts_with("study")),
            funs(case_when(is.na(.) ~ ".", 
                           TRUE ~ .))) %>%
  kable() %>% 
  kable_styling() %>%
  collapse_rows(1:3)

r1_spev_pv <- lm(spirit_score_std ~ por_score_std, d1)
r1_spev_abs <- lm(spirit_score_std ~ abs_score_std, d1)

r2_spev_abs <- lm(spev_score_std ~ abs_score_std, d2)
r2_dse_abs <- lm(dse_score_std ~ abs_score_std, d2)

r3_spev_por <- lm(spirit_score_std ~ por_score_std, d3)

r4_spev_por <- lm(spev_score_std ~ por_score_std, d4)
r4_dse_por <- lm(dse_score_std ~ por_score_std, d4)
r4_spev_pv <- lm(spev_score_std ~ pv_score_std, d4)
r4_dse_pv <- lm(dse_score_std ~ pv_score_std, d4)
r4_spev_abs <- lm(spev_score_std ~ abs_score_std, d4)
r4_dse_abs <- lm(dse_score_std ~ abs_score_std, d4)

df1b <- data.frame(study = "study 1",
                   outcome = "spiritual events",
                   predictor = c("porosity vignettes", "absorption"),
                   rsq = c(rsquared(r1_spev_pv)$R.squared, 
                           rsquared(r1_spev_abs)$R.squared))

df2b <- data.frame(study = "study 3",
                   outcome = c("spiritual events", "DSE"),
                   predictor = "absorption",
                   rsq = c(rsquared(r2_spev_abs)$R.squared, 
                           rsquared(r2_dse_abs)$R.squared))

df3b <- data.frame(study = "study 2",
                   outcome = "spiritual events",
                   predictor = "porosity scale",
                   rsq = c(rsquared(r3_spev_por)$R.squared))

df4b <- data.frame(study = "study 4",
                   outcome = rep(c("spiritual events", "DSE"), 3),
                   predictor = c(rep("porosity scale", 2), 
                                 rep("porosity vignettes", 2),
                                 rep("absorption", 2)),
                   rsq = c(rsquared(r4_spev_por)$R.squared, 
                           rsquared(r4_dse_por)$R.squared,
                           rsquared(r4_spev_pv)$R.squared, 
                           rsquared(r4_dse_pv)$R.squared, 
                           rsquared(r4_spev_abs)$R.squared,
                           rsquared(r4_dse_abs)$R.squared))

df_allb <- full_join(df1b, df2b) %>% full_join(df3b) %>% full_join(df4b) %>%
  mutate(outcome = factor(outcome,
                          levels = c("spiritual events", "DSE")),
         predictor = factor(predictor,
                            levels = c("porosity vignettes", "porosity scale", "absorption"))) %>%
  select(predictor, outcome, study, rsq) %>%
  arrange(predictor, outcome, study)

df_allb %>% 
  full_join(df_all %>% 
              filter(var == "spiritual experience") %>%
              select(-var) %>%
              rename(outcome = scale) %>%
              mutate(predictor = "country")) %>%
  mutate(predictor = factor(predictor,
                            levels = c("country", "porosity vignettes",
                                       "porosity scale", "absorption"))) %>%
  mutate(percent_exp = paste0(round(rsq * 100), "%")) %>%
  select(-rsq) %>%
  spread(study, percent_exp) %>%
  mutate_at(vars(starts_with("study")),
            funs(case_when(is.na(.) ~ ".", 
                           TRUE ~ .))) %>%
  kable() %>% 
  kable_styling() %>%
  collapse_rows(1:2)

Relationships

All studies, multipart plot

Spiritual Events only, cowplot

fig_s1_por <- d_all %>%
  filter(study == "Study 1") %>%
  ggplot(aes(x = pv_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Spiritual Events",
       color = "Country")
# fig_s1_por

fig_s1_abs <- d_all %>%
  filter(study == "Study 1") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s1_abs

fig_s1_title <- ggdraw() + 
  draw_label("STUDY 1", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s1 <- plot_grid(
  fig_s1_title,
  plot_grid(fig_s1_por, fig_s1_abs, ncol = 1, labels = c("A", "B")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s1

fig_s1_title_vert <- ggdraw() + 
  draw_label("STUDY 1", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s1_vert <- plot_grid(
  fig_s1_title_vert,
  plot_grid(fig_s1_por, fig_s1_abs, ncol = 2, labels = c("A", "B")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s1_vert

fig_s2_abs <- d_all %>%
  filter(study == "Study 2") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s2_abs

fig_s2_title <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s2 <- plot_grid(
  fig_s2_title,
  plot_grid(NULL, fig_s2_abs, ncol = 1, labels = c("", "C")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s2

fig_s2_title_vert <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s2_vert <- plot_grid(
  fig_s2_title_vert,
  plot_grid(NULL, fig_s2_abs, ncol = 2, labels = c("", "C")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s2_vert

fig_s3_por <- d_all %>%
  filter(study == "Study 3") %>%
  ggplot(aes(x = por_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Spiritual Events",
       color = "Country")
# fig_s3_por

fig_s3_title <- ggdraw() + 
  draw_label("STUDY 2", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 105))

fig_s3 <- plot_grid(
  fig_s3_title,
  plot_grid(fig_s3_por, NULL, ncol = 1, labels = c("D", "")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s3

fig_s3_title_vert <- ggdraw() + 
  draw_label("STUDY 2", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s3_vert <- plot_grid(
  fig_s3_title_vert,
  plot_grid(fig_s3_por, NULL, ncol = 2, labels = c("D", "")),
  ncol = 1, rel_heights = c(1, 10))
# fig_s3_vert

fig_s32_vert <- plot_grid(
  plot_grid(
    fig_s3_title_vert, 
    plot_grid(fig_s3_por, labels = c("C")), 
    ncol = 1, rel_heights = c(1, 10)),
  plot_grid(
    fig_s2_title_vert,
    plot_grid(fig_s2_abs, labels = c("D")), 
    ncol = 1, rel_heights = c(1, 10))
)

fig_s4_por1 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = pv_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_por1

fig_s4_por2 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = por_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_por2

fig_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = abs_score, y = spev_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Spiritual Events",
       color = "Country")
# fig_s4_abs

fig_s4_title <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0.5) +
  theme(plot.margin = margin(0, 0, 0, 190))

fig_s4 <- plot_grid(
  fig_s4_title,
  plot_grid(plot_grid(fig_s4_por1, fig_s4_por2, ncol = 2, labels = c("E", "F")), 
            plot_grid(NULL, fig_s4_abs, NULL, ncol = 3, rel_widths = c(1, 2, 1), labels = c("", "G", "")), 
            ncol = 1),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4

fig_s4_title_vert <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s4_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(plot_grid(fig_s4_por1, fig_s4_por2, ncol = 1, labels = c("E", "F")), 
            plot_grid(NULL, fig_s4_abs, NULL, ncol = 1, rel_heights = c(1, 2, 1), labels = c("", "G", "")), 
            ncol = 2),
  ncol = 1, rel_heights = c(1, 20))
# fig_s4_vert

fig_legend <- get_legend(fig_s4_por1 + theme(legend.position = "bottom"))

fig_all <- plot_grid(fig_s1, fig_s2, fig_s3, fig_s4, ncol = 4,
                     rel_widths = c(1, 1, 1, 2), scale = 0.95)
fig_all

fig_all_vert <- plot_grid(fig_s1_vert, fig_s32_vert, fig_s4_vert, fig_legend,
                          ncol = 1, rel_heights = c(1, 1, 2, 0.2))
fig_all_vert

Daily Spiritual Experiences only, cowplot

fig_s2_abs <- d_all %>%
  filter(study == "Study 2") %>%
  ggplot(aes(x = abs_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s2_abs

fig_s2_title_vert <- ggdraw() + 
  draw_label("STUDY 3", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s2_vert <- plot_grid(
  fig_s2_title_vert,
  plot_grid(fig_s2_abs, ncol = 1, labels = "B"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s2_vert

fig_s4_por1 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = pv_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Vignettes",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_por1

fig_s4_por2 <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = por_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Porosity Scale",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_por2

fig_s4_abs <- d_all %>%
  filter(study == "Study 4") %>%
  ggplot(aes(x = abs_score, y = dse_score)) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  xlim(0, 1) +
  ylim(0, 1) +
  theme(legend.position = "none") +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Absorption",
       y = "Daily Spiritual Experiences",
       color = "Country")
# fig_s4_abs

fig_s4_title_vert <- ggdraw() + 
  draw_label("STUDY 4", fontface = 'bold', x = 0, hjust = 0) +
  theme(plot.margin = margin(0, 0, 0, 7))

fig_s4_por1_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_por1, ncol = 1, labels = "A"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_por1_vert

fig_s4_por2_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_por2, ncol = 1, labels = "C"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_por2_vert

fig_s4_abs_vert <- plot_grid(
  fig_s4_title_vert,
  plot_grid(fig_s4_abs, ncol = 1, labels = "D"),
  ncol = 1, rel_heights = c(1, 10))
# fig_s4_abs_vert

fig_legend <- get_legend(fig_s4_por1 + theme(legend.position = "bottom"))

fig_all_vert <- plot_grid(
  plot_grid(fig_s4_por1_vert, fig_s2_vert, 
            fig_s4_por2_vert, fig_s4_abs_vert,
            ncol = 2),
  fig_legend,
  ncol = 1, rel_heights = c(2, 0.2))
fig_all_vert

Other versions

Spiritual Events only, one grid, new layout

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  gather(pred_scale, pred_score, c(por_score, pv_score, abs_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption"),
         pred_type = case_when(pred_scale == "Absorption" ~ "Absorption",
                               grepl("Porosity", pred_scale) ~ "Porosity",
                               TRUE ~ NA_character_),
         pred_type = factor(pred_type, levels = c("Porosity", "Absorption")),
         study_scale = paste(study, pred_scale, sep = ": "),
         study_scale2 = case_when(
           study == "Study 4"  & pred_scale != "Absorption" ~ 
             paste(study, pred_scale, sep = ": "),
           TRUE ~ study),
         study_scale3 = case_when(
           study == "Study 4" & pred_scale == "Porosity Scale" ~ "Porosity Scale",
           study %in% c("Study 1", "Study 4") ~ "Porosity Vignettes",
           study == "Study 2" ~ "Porosity Scale", 
           TRUE ~ " "),
         study_scale3 = factor(study_scale3,
                               levels = c("Porosity Vignettes",
                                          "Porosity Scale", " "))) %>%
  filter(!is.na(pred_score),
         spirit_scale == "Spiritual Events") %>%
  ggplot(aes(x = pred_score, y = spirit_score)) +
  facet_grid(rows = vars(pred_type), cols = vars(study, study_scale3)) +
  # facet_grid(pred_type ~ study_scale3) +
  geom_point(data = . %>% distinct(study, study_scale3, country, 
                                   pred_type, pred_scale, pred_score,
                                   spirit_scale, spirit_score),
             aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on predictor scale (Porosity Vignettes, Porosity Scale, or Absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Spiritual Events only, by study

fig_s1 <- d_all %>%
  filter(study == "Study 1") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s2 <- d_all %>%
  filter(study == "Study 2") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s3 <- d_all %>%
  filter(study == "Study 3") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

fig_s4 <- d_all %>%
  filter(study == "Study 4") %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Spiritual Events only, full grid

d_all %>%
  distinct() %>%
  gather(pred_scale, pred_score, c(pv_score, por_score, abs_score)) %>%
  mutate(pred_scale = recode_factor(pred_scale,
                                    "pv_score" = "Porosity Vignettes",
                                    "por_score" = "Porosity Scale",
                                    "abs_score" = "Absorption")) %>%
  # mutate(study_scale = paste(study, pred_scale, sep = ": ")) %>%
  filter(!is.na(pred_score)) %>%
  ggplot(aes(x = pred_score, y = spev_score)) +
  facet_grid(pred_scale ~ study) +
  geom_point(aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm",
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  # theme(legend.position = "bottom", 
  #       axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  # guides(color = guide_legend(override.aes = list(alpha = 1))) +
  guides(color = F) +
  labs(x = "Score on predictor scale\n(porosity or absorption; rescaled to 0-1)",
       y = "Score on Spiritual Events scale (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

Broken down by predictor type and spiritual scale

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  gather(poros_scale, poros_score, c(por_score, pv_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         poros_scale = recode_factor(poros_scale,
                                     "pv_score" = "Porosity Vignettes",
                                     "por_score" = "Porosity Scale"),
         study_scale = paste(study, poros_scale, sep = ": ")) %>%
  filter(!is.na(poros_score)) %>%
  ggplot(aes(x = poros_score, y = spirit_score)) +
  facet_grid(spirit_scale ~ study_scale) +
  geom_point(data = . %>% distinct(study, study_scale, country, 
                                   poros_scale, poros_score,
                                   spirit_scale, spirit_score),
             aes(color = country), alpha = 0.1) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on porosity measure (rescaled to 0-1)",
       y = "Score on spiritual experience measure (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")

d_all %>%
  gather(spirit_scale, spirit_score, c(spev_score, dse_score)) %>%
  mutate(spirit_scale = recode_factor(spirit_scale,
                                      "spev_score" = "Spiritual Events",
                                      "dse_score" = "Daily Spiritual Experiences"),
         study_scale = paste(study, "Absorption scale", sep = ": ")) %>%
  filter(!is.na(abs_score)) %>%
  ggplot(aes(x = abs_score, y = spirit_score)) +
  facet_grid(spirit_scale ~ study_scale) +
  geom_point(data = . %>% distinct(study, study_scale, country, 
                                   abs_score,
                                   spirit_scale, spirit_score),
             aes(color = country), alpha = 0.2) +
  geom_smooth(aes(color = country), method = "lm", 
              lty = 2, size = 0.7, alpha = 0, show.legend = F) +
  geom_smooth(method = "lm", color = "black", alpha = 0.7) +
  scale_color_brewer(palette = "Dark2") +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  guides(color = guide_legend(override.aes = list(alpha = 1))) +
  labs(x = "Score on absorption measure (rescaled to 0-1)",
       y = "Score on spiritual experience measure (rescaled to 0-1)",
       # caption = "Solid black lines correspond to to the overall trend, collapsing across countries",
       color = "Country")